Welcome

# These are the libraries we will use today. Make sure to install them in your console in case you have not done so previously.

library(tidyverse) # To load the collection of packages in the tidyverse
library(dplyr) # To wrangle our data (this package is also loaded by library(tidyverse))
library(readr) # To load the .csv data (this package is also loaded by library(tidyverse))
library(ggplot2) # To create plots (this package is also loaded by library(tidyverse))
library(tidyr) # To use pivot_wider() (this package is also loaded by library(tidyverse))
library(broom) # To format regression output
library(modelsummary) # To format model output
library(knitr) # To create HTML tables with kable()
library(kableExtra) # To format the HTML tables

library(janitor) # To examine our data with tabyl()
library(estimatr) # To run 2SLS with iv_robust()

Measuring the effect of mosquito net use on malaria infections

Imagine that the organization you work for is laying out a project to distribute mosquito nets to help combat malaria transmissions.

The funding agency requires a impact evaluation report from your organization. You are in charge of running the evaluation of this program.

You realize that there are potential unobserved confounders that could bias the observed differences in malaria risk for mosquito net users and non-users. You also think about the ethical considerations of fully randomizing who receives the nets, so you remember your Statistics II lecture on IVs and set up an encouragement design.

Your beneficiaries are scattered across ten villages. You decide to randomly select five villages to send SMS reminders every night encouraging them to use the mosquito nets. (This example is using simulated data)

Exploring our data

evaluation_df <- readr::read_csv("https://raw.githubusercontent.com/seramirezruiz/stats-ii-lab/master/Session%205/data/evaluation_data.csv") # loading simulated data frame of the intervention

You receive the results of your intervention from the M&E officers. There are 1000 inhabitants across the ten villages. This is what the data look like:

• village_name: A character string with the name of the village
• sms: A binary marker for the SMS encouragement
• net_use: A binary marker for mosquito net use
• malaria: A binary marker for malaria infection

Compliance types

You may remember, this table from the lecture:

We can cross-tabulate our data with janitor::tabyl() and the additional features of the janitor::adorn_* functions.

Why janitor::tabyl()? Because as prospective policy analysts we will do a lot of counting.

As the vignette of the package even puts it:

Analysts do a lot of counting. Indeed, it’s been said that ‘data science is mostly counting things.’ But the base R function for counting, table(), leaves much to be desired:

• It doesn’t accept data.frame inputs (and thus doesn’t play nicely with the %>% pipe)
• It doesn’t output data.frames
• Its results are hard to format. Compare the look and formatting choices of an R table to a Microsoft Excel PivotTable.

This is how it works. Say we are interested in exploring the number of persons in each of the observed strata, we would do:

evaluation_df %>% # your data frame
  janitor::tabyl(net_use, sms) %>% # the two dimensions for the table (D, Z)
  janitor::adorn_totals(c("row", "col")) %>% # add totals for rows and cols
  knitr::kable() %>% # turn into a kable table for nice rendering in HTML
  kableExtra::kable_styling() %>% 
  kableExtra::add_header_above(c("", "sms" = 2, "")) #add header for sms

If you want to learn more about the syntax of tabyl(), make sure to check the vignette

Exercise

Let’s explore the compliance types from this table

• Where are our compliers and non-compliers?
• How many people were encouraged via SMS, but did not use the net?
• How many people were not encouraged via SMS, yet they utilized the net?

Average malaria infections across strata

We can utilize the tabyl() syntax and our knowledge from the grammar of graphics to table and visualize the distribution of malaria on each stratum:

Table: Count of malaria infections across (Y,Z)

evaluation_df %>% # your data frame
  janitor::tabyl(malaria, sms) %>% # the two dimensions for the table (Y, Z)
  janitor::adorn_totals(c("row", "col")) %>% # add totals for rows and cols
  knitr::kable() %>% # turn into a kable table for nice rendering in HTML
  kableExtra::kable_styling() %>% 
  kableExtra::add_header_above(c("", "sms" = 2, "")) #add header for sms

Plot: Distribution of malaria infections across strata

ggplot(evaluation_df, aes(x = factor(sms), 
                          y = factor(net_use), 
                          color = factor(malaria))) +
  geom_jitter() +
  theme_minimal() +
  scale_x_discrete(labels = c("SMS = 0", "SMS = 1")) +
  scale_y_discrete(limits = c("1","0"), labels = c("NET = 1", "NET = 0")) +
  scale_color_manual(values = c("#CDCDCD", "#CC0055"),
                     labels = c("Not infected", "Infected")) +
  labs(x = "Encouragement",
       y = "Treatment",
       color = "") +
  theme(legend.position = "bottom")

Exercise

• What insights can we gather from the table and plot?

Exploring our set-up

lm(net_use ~ sms, data = evaluation_df) %>%
  modelsummary::modelsummary(stars = c('*' = .1, '**' = .05, '***' = .01),
                             statistic = 'conf.int',
                             fmt = 2)

naive_model <- lm(malaria ~ net_use, data = evaluation_df)

modelsummary::modelsummary(naive_model,
                           stars = c('*' = .1, '**' = .05, '***' = .01),
                           statistic = 'conf.int',
                           fmt = 2)

itt_model <- lm(malaria ~ sms, data = evaluation_df)
modelsummary::modelsummary(itt_model,
                           stars = c('*' = .1, '**' = .05, '***' = .01),
                           statistic = 'conf.int',
                           fmt = 2)

Let’s gather out local average treatment effect (LATE)

We have several options for this:

• Wald Estimator We can calculate this by hand. x. We can also calculate the average malaria rates among those who did and did not receive an SMS (no SMS = 0.504, yes SMS = 0.18).

Table: Observations across strata (D,Z)

Table: Count of malaria infections across strata (Y,Z)

Local Average Treatment Effect (LATE) manually

\[LATE = \frac{cov(Y_i,Z_i)}{cov(D_i,Z_i)}\]

Let’s take a look at our numerator \(cov(Y_i,Z_i)\), also \(ITT\)

\(E(Y|Z = 1) = \frac{90}{(410+90)} = 0.18\)
\(E(Y|Z = 0) = \frac{252}{(252+248)} = 0.504\)

Let’s take a look at our denominator \(cov(D_i,Z_i)\)

\(E(D∣Z = 1) = \frac{390}{(390 + 110)} = 0.78\)
\(E(D∣Z = 0) = \frac{165}{(165 + 335)} = 0.33\)

We can calculate our LATE

\[LATE = \frac{(0.18 - 0.504)}{(0.78 - 0.33)} = -0.72\]

• Two-stage Least Squares (2SLS). We will learn how to do this with iv_robust(), which is part of the estimatr package. It fits instrumental-variable regression through two-stage least squares. The syntax is the following:

estimatr::iv_robust(outcome ~ treatment | instrument, data)

In our case:

late_model <- estimatr::iv_robust(malaria ~ net_use | sms, data = evaluation_df)

modelsummary::modelsummary(late_model,
                           stars = c('*' = .1, '**' = .05, '***' = .01),
                           statistic = 'conf.int',
                           fmt = 2)

Exercise

• What is the substantive reading of these results?
• What would you tell the funding partner in your evaluation report?

Mechanics behind two-stage least squares (2SLS)

What estimatr::iv_robust() is doing in the background is the following:

net_use_hat <- lm(net_use ~ sms, evaluation_df)$fitted.values # get fitted values from first stage (the part of x that is exogenously driven by z)
summary(lm(evaluation_df$malaria ~ net_use_hat)) # run second stage with instrumented x (careful, the standard errors are wrong; better use iv_robust() from estimatr instead)

## 
## Call:
## lm(formula = evaluation_df$malaria ~ net_use_hat)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -0.504 -0.180 -0.180  0.496  0.820 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.74160    0.03757   19.74   <2e-16 ***
## net_use_hat -0.72000    0.06273  -11.48   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4463 on 998 degrees of freedom
## Multiple R-squared:  0.1166, Adjusted R-squared:  0.1157 
## F-statistic: 131.8 on 1 and 998 DF,  p-value: < 2.2e-16

Thinking about the validity of instruments

We can also adapt the iv_robust() syntax to accommodate valid conditional instruments:

estimatr::iv_robust(Y ~ D + W | Z + W, data = your_df)